#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <sys/time.h>
#include <tbb/tbb.h>

using namespace std;

typedef double FD;

template <typename T, typename F> T integrate(T st, T en, int div, F func) {
  T localRes = 0;
  T step = (en - st) / div;
  T x;
  x = st;
  localRes = func(st) + func(en);
  localRes /= 2;
  for (int i = 1; i < div; i++) {
    x += step;
    localRes += func(x);
  }
  localRes *= step;

  return localRes;
}

// template <typename T, typename F> T integrate_mp(T st, T en, int div, F func) {
//   T localRes = 0;
//   T step = (en - st) / div;
//   T x;
//   x = st;
//   localRes = func(st) + func(en);
//   localRes /= 2;

// #pragma omp parallel for private(x)
//   for (int i = 1; i < div; i++) {
//     x = st + i * step;
//     T temp = func(x);
// #pragma omp critical
//     localRes += temp;
//   }

//   localRes *= step;

//   return localRes;
// }

template <typename T, typename F>
T integrate_mp_reduce(T st, T en, int div, F func) {
  T localRes = 0;
  T step = (en - st) / div;
  T x;
  x = st;
  localRes = func(st) + func(en);
  localRes /= 2;

#pragma omp parallel for private(x) reduction(+ : localRes)
  for (int i = 1; i < div; i++) {
    x = st + i * step;
    localRes += func(x);
  }

  localRes *= step;

  return localRes;
}

template <typename T, typename F> T integrate_tbb(T st, T en, int div, F func) {
  T localRes = 0;
  T step = (en - st) / div;
  T x;
  x = st;
  localRes = func(st) + func(en);
  localRes /= 2;
  localRes += tbb::parallel_reduce(
      tbb::blocked_range<int>(0, div), localRes,
      [=](tbb::blocked_range<int> &r, T res) {
        for (int i = r.begin(); i < r.end(); ++i) {
          res += func(st + i * step);
        }
        return res;
      },
      std::plus<T>());

  localRes *= step;

  return localRes;
}

template <typename T, typename F>
T test(T st, T en, int div, F func, T (*prec)(T, T, int, F)) {
  struct timeval t_start, t_end;
  long timeuse = 0;
  gettimeofday(&t_start, NULL);

  FD Res = prec(st, en, div, func);

  gettimeofday(&t_end, NULL);
  timeuse = 1000000 * (t_end.tv_sec - t_start.tv_sec) + t_end.tv_usec -
            t_start.tv_usec;
  printf("time=%f\n", timeuse / 1000000.0);

  return Res;
}

int main(int argc, char *argv[]) {
  FD start = 1;
  FD end = 5;
  int div = INT32_MAX;
  auto func = [](FD x) { return x * x + 2 * sin(x); };

  FD Res_Seq = test(start, end, div, func, integrate);
  // FD Res_MP = test(start, end, div, func, integrate_mp);
  FD Res_MP_reduce = test(start, end, div, func, integrate_mp_reduce);
  FD Res_TBB = test(start, end, div, func, integrate_tbb);

  printf("Res_Seq   = %lf\n", Res_Seq);
  // printf("Res_MP    = %lf\n", Res_MP);
  printf("Res_MP_r  = %lf\n", Res_MP_reduce);
  printf("Res_TBB   = %lf\n", Res_TBB);
  return 0;
}